### this code was developed for the Eggers & Vivyan project. see pivotprobs package for updated version.

# simulation based
pivotal.probabilities = function(samp, tol = .01, out.as.vector){
	rank.mat = t(apply(-samp, 1, rank)) # why is there a transpose?
	out = matrix(NA, ncol(samp), ncol(samp))
	for(i in 1:ncol(samp)){
		for(j in 1:ncol(samp)){
			if(i >= j){next} # this will be a lot faster.
 			out[i,j] = sum( (rank.mat[,i] + rank.mat[,j]) == 3 & abs(samp[,i] - samp[,j]) < tol)
		}
	}
	out.mat = out/nrow(samp)
	if(!out.as.vector){return(out.mat)}
	data = as.vector(out.mat)
	labels = c(); for(name in colnames(samp)){labels = c(labels, paste0("p.", name, ".", colnames(samp)))}
	out = as.vector(na.omit(data))
	names(out) = labels[!is.na(data)]
	out	
}

win.shares = function(samp){
	max.of.row = apply(samp, 1, max, na.rm = T)
	out = c()
	for(j in 1:ncol(samp)){
		out = c(out, sum(samp[,j] == max.of.row))
	}
	names(out) = colnames(samp)
	out/sum(out)
}


### new analytical (naive) version

library(gtools)
# require(ExtDist)

# first note that pBeta_ab gives the same thing as pbeta when we rescale the variable. 
# y = .27
# alpha.vec = c(40, 30, 20, 10, 5)
# i = 3
# pBeta_ab(y, alpha.vec[i], sum(alpha.vec[-c(1,2,i)]), a = 0, b = 1 - 2*y)
# pbeta(y/(1 - 2*y), alpha.vec[i], sum(alpha.vec[-c(1,2,i)]))

napotffay = function(alpha.vec, y){
	# naive.analytical.probability.of.tie.for.first.at.y 
	# probability that parties 1 and 2 tie at y can be calculated without specifying the values of the other parties, by the aggregation property.
	prob.of.tie.at.y = ddirichlet(x = c(y,y,1-2*y), alpha = c(alpha.vec[1], alpha.vec[2], sum(alpha.vec) - alpha.vec[1] - alpha.vec[2]))
	# now, what is the probability that _no other party_ would have gotten above y?
	# we do it separately for each party -- a naive approach not recognizing how the vote shares for the other parties are related.
	# obvious they are correlated because they must sum to 1-2*y.  
	probs = c()
	for(i in 3:length(alpha.vec)){
	  this.prob = pbeta(y/(1-2*y), alpha.vec[i], sum(alpha.vec[-c(1,2,i)])) # this is the part I'm least sure of
# 		this.prob = pBeta_ab(y, alpha.vec[i], sum(alpha.vec[-c(1,2,i)]), a = 0, b = 1 - 2*y) # this is the part I'm least sure of
		probs = c(probs, this.prob)
	}
	list(prob.of.tie.at.y = prob.of.tie.at.y, probs.of.being.below.y = probs, probability = prod(c(prob.of.tie.at.y, probs))) 
}

# now we aggregate this up (integrate over values of y)
probability.of.tie.for.first = function(alpha.vec, increments = 50){
	ys = seq(1/length(alpha.vec), .5, length = increments) # the least a pair of parties can get and be tied for first is 1/k each. the most they can get is .5.  We are going to calculate the probability of a tie for first at each value in that sequence and then integrate.
	probs = c() # this holds the probability of the two parties being tied for first at each y
	tie.vec = c()  # this holds the probability of the two parties being tied at a given y 
	prob.mat = matrix(NA, nrow = 0, ncol = length(alpha.vec)-2) # this holds the probabilities of being below y for each of the parties other than the ones who are tied.
	increment = ys[2] - ys[1]
	for(i in 1:(length(ys) - 1)){
		this = napotffay(alpha.vec, (ys[i] + ys[i+1])/2)
		probs = c(probs, increment*this$probability) 
		tie.vec = increment*c(tie.vec, this$prob.of.tie.at.y)  # this seems wrong, though it's not used. 
		prob.mat = rbind(prob.mat, this$probs.of.being.below.y) 
	}
	list(estimate = sum(probs), # apply(cbind(tie.vec, prob.mat), 1, prod)),  # for each y, the probability of 
		ys = ys, tie.vec = tie.vec, prob.mat = prob.mat) # I include these for diagnostic purposes
}

## and here's a function that yields pivotal probabilities in a vector, given an alpha input.
pivotal.probabilities.analytical = function(alpha.vec, increments = 50, normalize = F, as.matrix = F, n = 50000){
	out.vec = name.vec = c()
	out.matrix = matrix(0, nrow = length(alpha.vec), ncol = length(alpha.vec))
	rownames(out.matrix) = colnames(out.matrix) = names(alpha.vec)
	for(i in 1:length(alpha.vec)){
		for(j in 1:length(alpha.vec)){
			if(i <= j){next}
			this.entry = probability.of.tie.for.first(alpha.vec = c(alpha.vec[i], alpha.vec[j], alpha.vec[-c(i,j)]), increments = increments)$estimate/n
			out.matrix[i,j] = out.matrix[j,i] = this.entry
			out.vec = c(out.vec, this.entry)
			name.vec = c(name.vec, paste0("p.", names(alpha.vec)[i], ".", names(alpha.vec)[j]))
		}
	}
	if(normalize){
		out.vec = out.vec/max(out.vec)		
	}
	names(out.vec) = name.vec
	if(as.matrix){return(out.matrix)}
	out.vec
}



